Stochastic simulations of fermionic dynamics with phase-space representations 



O 

(N 



M. Ogren, K. V. Kheruntsyan, J. F. Corney 
ARC Centre of Excellence for Quantum-Atom Optics, School of Mathematics and Physics, 
University of Queensland, Brisbane, Queensland 4072, Australia 
(Dated: March 17, 2011) 

A Gaussian operator basis provides a means to formulate phase-space simulations of the real- and 
imaginary-time evolution of quantum systems. Such simulations are guaranteed to be exact while the 
underlying distribution remains well-bounded, which defines a useful simulation time. We analyse 
the application of the Gaussian phase-space representation to the dynamics of the dissociation of an 
ultra-cold molecular gas. We show how the choice of mapping to stochastic differential equations can 
be used to tailor the stochastic behaviour, and thus the useful simulation time. In the phase-space 
approach, it is only averages of stochastic trajectories that have a direct physical meaning. Whether 
particular constants of the motion are satisfied by individual trajectories depends on the choice of 
mapping, as we show in examples. 
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Numerical approaches are an indispensable part of en- 
deavours to understand quantum many-body physics in 
condensed matter and AMO physics. In particular, there 
is a need for real-time, dynamical simulations, driven in 
large part by the progress in the control and flexibility of 
ultra-cold atom experiments, which has made the dynam- 
ically evolving quantum many-body state more directly 
accessible. For bosons, first-principles phase-space meth- 
ods have successfully simulated dynamics in experimen- 
tally realistic systems [3, HI- However, these methods are 
not directly applicable to fermionic systems, which are 
an increasingly important area of ultra-cold atoms, often 
with direct relevance to condensed matter systems. 

The exponential growth of the Hilbert space with sys- 
tem size hinders a brute-force approach for systems of 
more than a few modes. Stochastic approaches, provided 
they are unbiased, can provide exact results within the 
precision determined by sampling error. A range of quan- 
tum Monte Carlo methods has been used to address a va- 
riety of problems in many-body quantum physics. How- 
ever, the limitations when it comes to dynamics are well 
known Q, for example, the oscillating phase problem in 
path- integral approaches An interesting direction in 
recent years has been the extension to fermionic systems 
of stochastic wavefunction approaches 

In this work we employ a Gaussian stochastic method 
based on a generalized phase-space representation of the 
quantum density operator Q . The representation allows 
the quantum Liouville equation for the density operator 
to be mapped onto an equivalent Fokker-Planck equa- 
tion for a distribution function over phase space, provided 
that the distribution vanishes at the boundary. This dis- 
tribution is then sampled via equivalent stochastic phase- 
space equations, with physical results corresponding to 



stochastic averages. The phase-space equations are struc- 
turally similar to the Heisenberg equations for the corre- 
sponding operators, with additional stochastic noise. 

Here we explore the freedom in choosing the stochas- 
tic noise in order to reduce sampling errors and extend 
the useful simulation time. We first introduce the Bose- 
Fermi model we use study these effects, and a set of con- 
served quantities that can be used to benchmark the dif- 
ferent choices of stochastic equations. After reviewing the 
phase-space formalism, we give the general form of the 
stochastic equations corresponding to the Hamiltonian. 
To exemplify the gauge freedom, we then give two forms 
of the noise terms and demonstrate through simulation 
their very different numerical properties. 

As a particular application, we consider a model of 
production of correlated pairs of fermionic atoms by dis- 
sociation of a Bose-Einstein condensate (BEC) of di- 
atomic molecules @, 3- The uniform molecular BEC 
is initially in a coherent state at zero temperature with 
average initial number of molecules Nq, with no atoms 
present. The created fermionic atoms are modeled as be- 
ing untrapped in the x-direction, and propagate through 
the homogeneous condensate. The Hamiltonian of this 
boson-fermion model [9] is given by 



H = hy Ak^k a — ihn 



(1) 



where k labels the M plane-wave modes for a quantiza- 
tion box of length L and a — 1 , 2 labels the effective spin 
state for the atoms. The fermionic number and pair oper- 
ators are defined as nk,o- = c k a c\. a and = Ck,ic_k.2, 
respectively, with {ck, CT , cj^, a , } — #k,k'£a-,o-'j while the 
bosonic molecular operators obey [a, a*] = 1. The first 
term gives the kinetic energy of the atoms of mass m a 
and the detuning A between the atomic and molecular 
levels: hA^ = ti 2 |k| 2 /(2m a ) + HA. The second term 
describes the atom-molecule coupling of strength k . 

The physics of the growth of correlations during the 
dynamics has been explored elsewhere Here, we use 
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the same system parameters but focus on the evolution of 
certain conserved quantities. While such quantities are 
constant in the stochastic averages, which have a phys- 
ical meaning, they are not necessarily constant in indi- 
vidual trajectories. We study this evolution to monitor 
the growth of sampling error for different choices of the 
stochastic equations, and to illustrate the exactness of 
the method within the limitations of sampling error. 

The spin-symmetry of the Hamiltonian implies the 
identity = n.k,i = n._ k,i = "-k,2 = "-k,2 for equal ini- 
tial populations. An additional operator identity arises 
from the homogeneity of the molecular condensate, 



m^TOk (= n k ,i"--k,2) = "-k- 
According to this, we expect the conserved quantity 

F k = (m\jh-kj - (n k ) 



(2) 



(3) 



to be zero in any numerical implementation. We calcu- 
late this quantity numerically for the resonant Fourier 
mode ko, along with the total energy normalised by the 
dissociation energy: E = (H)/2h\A\ and the total num- 
ber of molecules and pairs, normalised by the initial num- 
ber molecules: N = (2(a^a) + J2w afok^)) /2No, which 
is also conserved. 

The Gaussian phase-space representation maps pairs 
of annihilation/creation operators onto first-order differ- 
ential operators. It can thereby be used to transform the 
Liouville equation for unitary evolution 



dt P h 



(4) 



into a differential equation for an equivalent phase-space 
distribution, so long as certain boundary terms vanish. In 
practice the appearance of a boundary term is indicated 
by the rapid growth of sampling error and the appearance 
of large excursions in the trajectories [Hj|, an( i this places 
a limitation on the length of the simulation. For Hamil- 
tonians containing up to four operators, a second-order 
partial differential equation is generated, which can be 
written in the form of a Fokker-Planck equation (FPE) : 



dt 



P(A) 



8 



i J 



1 



d 2 



2 ^— dXjdXk 

j,k J 



Dij(X) 



P(A). 
(5) 

The first order derivatives in the phase-space variables 
Xj correspond to drift behaviour and the second order to 
the diffusion. Effects such as three-body interactions will 
result in higher-order derivatives, but these are difficult 
to efficiently sample by numerical methods (l2^ . 

In general the phase-space variables A are complex. 
However, the analytic nature of the Gaussian operators 
gives a freedom in the choice of derivatives when the vari- 
ables are expanded into real and imaginary parts. The 



diffusion matrix D of the resulting FPE can always be 
chosen to be positive-definite [13], as required for stochas- 
tic sampling. 

The quantum state generated by the Hamiltonian ([T]) 
can be represented by a distribution over 3M+2 variables 
X(t) = (m, . . . ,HM,mi, • • ■ ,m A /,mf , . . . , , /3, /3+), 
with =^ m* and f3 + ^ (3*. The corresponding FPE 
for the dynamics is 

d t P = 2i J2k A k <9 mk m k - d m + m i p 

+«Ek (P +m * + PK) - d ™ k P (1 - 2n k ) 

-d m +/ 3+ ( 1 - 2n k) + c^m k + d,3 + TO+ (6) 



d m + d P n l + d m +dp + m 



Note that all differential operators act also on the multi- 
dimensional distribution P = P (S>, tj . To directly solve 

Eq. ([6]) is computationally unfeasible for many variables. 
Instead one can employ a mapping [lil [l5| to an equiv- 
alent set of stochastic differential equations (SDEs) to 
sample the moments of the distribution. In the Ito cal- 
culus, stochastic equations corresponding to Eq. ([5]) have 
the general form 



dnu = (( 



"m k ) dr + A(7 1/2 B( nk )dW 



drrik = [— 2i£k?7!.k + < 
drri^ = [2if>k^k + Q 

da = _ ]^o Ek TO k dr 

A* + = -^Ek 



(1- 
(1- 



2n k )] dr 4 
2nk)l dr ■ 



-1/2 



B^dW, 
1/2 BW)dW, 



-1/2 



B^dW, 



midr 



l/2 B ( Q - t 



>dW, 



(7) 

where we have used a scaled time, r = KyNot and have 
also normalized the molecular field by its maximum (ini- 
tial) value, i.e. a = /3/y/ Nq . The deterministic part 
of the Ito equations corresponds to the drift terms in the 
FPE, which if taken alone, are equivalent to the so-called 
'pairing mean-field theory' 0, 0, E3|- The stochastic 
part, in which B"' are row vectors with elements that are 
functions of the phase-space variables, and where dW is a 
column- vector of real Wiener increments, constitutes dif- 
fusion processes in the complex phase-space. This form 
of Eqs. ([7]) shows that with drift terms of order 1, the 
stochastic terms are of order l/ v / Ao, i.e. the stochastic 
terms and therefore non-mean-field corrections are more 
important for decreasing No. 

Stochastically sampled moments can be related to 
physical expectation values. For example, the first or- 
der moments give: 




(m k ) = (&k,i&-k,2) 



(8) 



(a)/VNo- 



Normally ordered higher-order moments are obtained ex- 
actly by stochastic averages of a corresponding Wick de- 
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composition as in the following example 



(m+m k ) s + (ra£)s = (m k ra k ) 



(9) 



Note, however, that this does not mean that a Wick 
factorisation is assumed to hold for a general quantum 
state [l0( , since the average of a product is not the same 
as the product of averages. 

The equivalences above hold so long as the appropriate 
moments of the distribution are well-defined. In prac- 
tice this requires that the distribution tails vanish suffi- 
ciently quickly, which again places a limit on the simula- 
tion time, indicated by 'spiking' behaviour and associated 
rapid growth of sampling error. The instabilities under- 
lying this behaviour are a general feature of nonlinear 
stochastic equations (ll. 18|. 

With the equivalences between stochastic averages and 
operator expectation values, the defined conserved quan- 
tities Fk [Eq. |3])], E and N can be calculated: 



E 



F k = 



k"k 



(m^TOk + rt k - n k )s, 



(10) 
(11) 

(12) 



Note that although the stochastic quantities defined in 
Eqs. (fTU|) - (fT2l are complex for individual trajectories, the 
average of the imaginary components approach zero with 
increasingly many stochastic trajectories sampled. Thus 
the average of each of these quantities approaches a real 
value, as expected for physical observables. 

The stochastic equations corresponding to a given 
Hamiltonian are not unique and therefore can be tai- 
lored to give different numerical and sampling proper- 
ties 14|. We illustrate how this can be done through 
the choice of 'diffusion gauges' to extend the useful sim- 
ulation time 0, [3, [2(| . The stochastic terms must ful- 
fill the matrix-square-root condition [3] that relates the 
diffusion matrix D in the Fokker-Planck equation to the 
noise-matrix B : 



D = BB , B = 



B("k) B (mk) . B( m k) B (a) b( q+ ) 



(13) 

where T denotes matrix transpose. Let O denote a ma- 
trix with orthonormal rows composed of functions of 
phase-space variables. Then if B fulfills Eq. (fT3f . so does 
B = BO, which gives infinitely many choices of the SDE. 

One specific noise matrix, which we together with 
Eq. label SDE-I, is 



Br = 



ml 
-nl 


—iml 
inl 
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(14) 
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Figure 1: Real and imaginary part of the conserved quan- 
tity Jk , defined in Eq. (I10|l . as a function of scaled time t. 
Dashed line gives the mean for SDE-I, while the sampling er- 
ror vanishes for SDE-I (see text). Solid line gives the mean for 
SDE-II, with the light-grey shading giving the sampling error 
(±<t). The spiking behaviour and rapid growth of sampling 
error at r = 1.5 mean that the results from SDE-II cannot be 
used past this time. The insets show the standard deviation 
a of Fk g for SDE-II, with the arrows pointing out precursors 
of the spiking behaviour (see text). We use a momentum grid 
with |fc] = 0, dk, lOOOdfc, dk = 2-k/L ~ 7.92 x 10 3 m _1 
and a resonance momenta ko = ^j2m a \A\/h = 500dk, where 
A = —12500 s . Initially we have No = 100 molecules and 
an atomic vacuum. The atom-molecule coupling strength is 
k ~ 500 s - . The stochastic quantities for both SDE-I and 
SDE-II are evaluated using 10 4 trajectories. 



where dWj = [ dw\ dw2 dw^ dw4 ] T 
Note that it is often notationally convenient 
to work instead with complex Wiener incre- 
ments, e.g. dW {1) = {dwi +idw 2 )/V2, such 
that dW (]) satisfies (dW^ (t) dwO') (t')) 
0, (dW^ (t) dw( j ')* (t')) = S jf S (r - r') dr. Then, for 
example, B^ k) dW/ = n k (m^dW {l > + m+dW {2 >) . 

As proved in Appendix A, this choice of noise terms 
means that the quantity Fk defined in Eq. ([TO")) is satisfied 
by each individual trajectory, not just by the ensemble 
average, i.e. 



mjm k 



n k . 



(15) 



This property is clearly seen graphically in Fig. fT] as a 
vanishing sampling error for SDE-I. In Figs. [2] and [31 we 
see that the energy and particle number are conserved 
for SDE-I, but with a finite sampling error (dark-grey 
shading) which can be reduced further by including more 
stochastic trajectories. The trajectories are stable, with 
no 'spiking' or dramatic increase in sampling error, un- 
til at least a normalised time of t = 5.0. We conclude 
that SDE-I performs very well for the particular set of 
parameters chosen. 

We now use the gauge freedom of the condition in 
Eq. ([13")) to construct another specific noise matrix (SDE- 
II) which does not fulfil Eq. (|15[) . For the case of a single 
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Figure 2: Normalised total energy E, denned in Eq. 
as a function of scaled time r. Dashed line and dark-grey 
shading give the mean and sampling error (±cr) for SDE-I. 
Solid line and light-grey shading give the mean and sampling 
error (±<r) for SDE-II. While the sampling error grows in both 
cases, SDE-I is stable for at least 3 times longer. Parameters 
are as in Fig. [T] 



k-mode, the noise matrix for this diffusion gauge can be 
written: 







-ink 
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(16) 

However, in general Bjj is of size (3M + 2) x 8M, 
i.e. the number of noise columns grows with the num- 
ber of phase-space variables, such that now dW// = 

iT 



[ dwi^k dw 



2.k 



for example B 



dw 8t k J J / y/2 
'dW// = 



In this case we have 



»k 



Ek TO k ( 



dwP* + dW, 



(2)- 



and 



dW, 



(i) 



dW, 



(3) 



For SDE-II it is now only the average of that is 
zero, within the finite sampling error indicated by the 
light-grey shading in Fig. [TJ For the energy and particle 
number, the average is still constant within the sampling 
error, as shown in Figs. [5] and [3] However, the sampling 
error is now larger than for SDE-I. Moreover, the mean 
results from SDE-II (solid line in all figures) start to spike 
before r ~ 1.5, with an associated dramatic increase in 
sampling error, and can thus not be used beyond this 
point for the present parameters. The standard deviation 
of a stochastic variable is a moment of higher order than 
the average of the variable itself, and precursors of the 
spiking behaviour are first seen here. This is illustrated 
by the arrows in the inset plots of Fig. Q]for the standard 
deviation of i*k > but generally occurs for all sampled 
variables. 

As shown in Figs.[T][3l there can be dramatic differences 
in the performance of the two diffusion gauges. However 
the relative performance depends on the system param- 



Figure 3: Normalised total particle number as defined in 
Eq. (|12p . as a function of scaled time r. Dashed line and 
dark-grey shading give the mean and sampling error (±cr) for 
SDE-I. Solid line and light-grey shading give the mean and 
sampling error (±a) for SDE-II. Parameters are as in Fig. [T] 



eters. For instance, whereas the second gauge (SDE-II) 
may seem unnecessarily complicated for the many modes, 
and leads to a much larger sampling error and a shorter 
useful simulation time here, for the small system in [Io| . 
it is in fact superior to SDE-I in terms of useful simula- 
tion time. 

From the theoretical foundation it is expected that the 
Gaussian phase-space method is exact while the distri- 
bution is sufficiently bounded [6j. In practice the sim- 
ulation can be trusted until signatures such as spiking 
trajectories and rapid growth of the sampling error oc- 
curs in the time evolution of the phase-space variables 
[ToL [ill I20L HH ■ We have previously also analysed a re- 
lated dynamical system with only jVo = 10 molecules 
and M — 10 atomic momentum modes [10]. For this 
test system, the exponentially growing dimension of the 



Hilbert space was small enough (d = 2 



10 5 ), to 



allow a direct comparison to an expansion in a number 
state basis. However, this comparison is not possible for 
the system under study here. Having explicit access to 
different stochastic realisations of the FPE, as here with 
Eqs. (Q3J) and ([To]) , then gives the possibility to compare 
different stochastic calculations of the moments to check 
the accuracy of the numerical implementation or to de- 
tect errors in the underlying derivations. 

Despite the different stochastic behaviour revealed in 
Figs. OB it is important to note that SDE-I and SDE- 
II both correspond to the same Hamiltonian ([lj and the 
same complex FPE Eq. ((6j). Underlying these different 
realisations is the overcompleteness of the Gaussian rep- 
resentation, which allows the one density operator p to 
be mapped to many different distributions. 

In summary, we have demonstrated how different diffu- 
sion gauges can substantially change the numerical per- 
formance of the Gaussian fermionic phase-space method. 
This ability to manipulate the form of stochastic equation 
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can be used to reduce the sampling error and extend the 
useful simulation time, depending on the system param- 
eters. In addition, we have shown that the simulation 
of conserved quantities can have qualitatively different 
behaviour for different gauges. The conserved quantities 
thus provide a check on numerical implementation and 
allow the performance of different gauges to be bench- 
marked. 
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APPENDIX A. DERIVATION OF EQ. (Tig]) 

Here we prove Eq. (fl5| for SDE-I, which is a stronger 
condition than the corresponding result for the stochastic 
average. We apply the product rule for two stochastic 
variables X and Y within Ito calculus 

di (XY) = Xdi (Y) + di (X) Y + di (X) di (Y) , (17) 

to the first term in Eq. (|15|l . with di denoting the Ito 
differential. Hence we have, from Eqs. (0 and (|T4|) 



dj (mjmk) = — 2i5\ <i m^m^dT + am^ (1 — 2nk) dr 
(mldWl - n k dW 2 *) + 2i5^m+m^dT 



+N, 



+a+ (1 - 2n k ) m k dr + N 1/2 (m+ 2 dW^ - n^dW?) m k 

(m£ 2 dW£ - nldWl) (m£dW? - n\dW^) 
— {am^ + a + rnk) (1 — 2nk) dr 



+N f 



-1/2 



:m k - nQ (m k dW* + m£dW%) , 



(18) 

where we have kept, as usual, only first order terms in 
dr. The increment for n 2 , can be calculated similarly, 
leading to the following equation for the increment in 
the left-hand side of Eq. f[TB]>: 



di (m£ TOk + = 
+N- 1/2 (m.' 



T-k^k 



(am^ + a + TOk) dr 
£) (m k dW? + m+dWi 



(19) 



From Eqs. (|7|) and ()14[) . the corresponding expression for 
the left-hand side of Eq. (JX5J is 

drik — (am^ + a + mk) dr+N Q 1 ^ 2 nk (m^dW* + m£dW£ ) 

(20) 
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